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I .   INTRODUCTION 

This  paper  brings  together  the  various  scattered  results  on  the 
Renewal  Process  model  which  are  most  relevant  to  applications.  The 
early  sections  are  expository  in  nature,  describing  the  model  and  its 
applications,  with  particular  emphasis  on  the  Poisson  process.  Sub- 
sequent sections  are  devoted  to  general  bounds  and  approximations  for  the 
renewal  intensity,  renewal  function  and  forward  recurrence  distribution. 
A  few  special  cases  are  worked  out  exactly  and  a  numerical  method  for 
the  discrete  case  is  given.  Special  approximations  for  the  Weibull 
case  are  included,  and  bounds  on  higher  moments  of  the  renewal  counting 
process  are  given.  Numerical  illustrations  are  given  whenever  helpful. 
A  theoretical  development  of  the  Renewal  Process  model  is  found  in 
[Smith,  1958].  The  entire  subject  matter,  including  some  extentions,  is 
fully  discussed  there  and  an  extensive  bibliography  is  provided.  A  good 
introduction  to  the  subject  is  [Cox,  1962],  where  the  fundamentals  are 
treated  in  a  \/ery   clear  manner,  for  the  case  where  inter-event  times 
have  a  density.  An  excellent  discussion  of  Renewal  Theory  is  in  [Feller,  1971], 
where  a  section  is  devoted  to  the  two-sided  case  as  well.  We  will  not  consider 
the  more  general  theory  based  on  two-sided  distributions,  as  our  applications 
call  for  the  nonnegative  case  only.  Some  references  to  this  generality  are 
[Stone,  1965  and  1966].  Proofs  of  most  results  are  not  given  here,  but  in- 
stead intuitive  reasoning  is  provided  to  assist  the  reader's  insight  and 
comprehension.  Finally,  the  problems  of  estimation  of  rates,  tests  for 
trends,  and  other  issues  of  a  statistical  nature  are  left  to  [Cox  and  Lewis, 
1966]  for  resolution.  This  book  gives  an  excellent  discussion  of  the  statistical 
problems  which  arise  in  practice. 


II.  POINT  PROCESS 

A  point  process  is  a  series  of  epochs,  t-,  ■  t?   t-  ...  in  time, 
which  describe  the  evolution  of  a  process,  usually  stochastic  in  nature. 
This  model  is  used  to  describe  such  processes  as  fail  ure  epochs  of  complex 
equipments,  emissions  of  particles  during  decay  or  pulses  along  a  nerve 
fibre  for  example.  The  renewal  process  is  a  special  kind  of  point  process. 
If  the  interevent  times  t,  ,  t?  -  t,  ,  t,  -   t? ,  . . .  ,  are  statistically 
independent  random  variables,  and  all  except  possibly  the  first  have 
a  common  probability  distribution,  we  have  a  renewal  process. 

The  simplifying  assumptions  of  a  renewal  process  make  it  a  frequent 
choice  as  a  model  for  certain  point  processes.  One  renewal  process  is  most 
often  assumed  and  its  properties  are  particularly  worth  discussing,  namely 
the  Poisson  process.   It  is  the  standard  or  benchmark  to  which  other  point 
processes  are  compared. 

III.  THE  POISSON  PROCESS 

Consider  a  point  process  whose  events  occur  at  the  epochs  denoted  by 

o  <  t}   <  t2 .<  t3  ...  .  if  x1  ?  t-,,  xz  =  tz  -  t1  ,  x3  =  t3  -  u  ,  ... 

are  statistically  independent  random  variables  each  with    the  exponential 
distribution, 

Prob(Xn  s  x}   =   1    -  e"Ax,   x  >_  0    , 

then  we  call  the  process  a  (homogenous)  Poisson  process.  The  X   are  called 

n 

interevent  times  and  the  parameter  A  ,  the  reciprocal  of  the  mean  (average) 
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time  between  events,  is  called  the  rate  of  the  process,  as  it  is  the  mean 
number  of  events  occurring  per  unit  time  over  any  time  interval.   In  fact, 
the  actual  number  of  events  occurring  per  unit  time,  over  an  increasing  time 
period,  tends  almost  surely  to  A  . 

There  are  many  physical  processes  for  which  evidence  suggests  the 
Poisson  process  as  an  appropriate  model.  These  include  telephone  calls 
made  during  a  particular  hour  of  the  day,  emissions  from  a  radioactive 
source  over  a  period  during  which  the  strength  of  the  source  remains  substan- 
tially constant,  the  occurrence  of  nerve  pulses  along  a  nerve  fibre,  or 
the  failures  of  complex  equipment  during  its  operating  periods.  See  in 
particular  [Drenick,  1960]  for  some  general  conditions  under  which  the 
failures  of  complex  equipment  form  approximately  a  Poisson  process.  Also, 
[Cox  &  Lewis,  1966,  p.  251]  show  some  data  from  several  physical  processes 
which  exhibit  the  Poisson  process  assumptions.  The  following  description 
is  particularly  suggestive  of  its  wide  applicability. 

Consider  an  arbitrary  but  fixed  time  horizon  T  ,  and  a  population 
of  size  n  of  individuals  who  may  place  a  call  during  (0,  T]  .  Suppose 
further  that  each  individual  decides  to  make  a  call  with  probability  p  , 
and  that  all  individuals  act  independently,  so  that  the  number  calling  in 
(0,  T]  has  the  Binomial  distribution  with  parameters  n,  p  .  Assume 
further  that  those  calling  select  a  time  to  call  independently  of  others 
and  uniformly  from  (0,  T]  .  The  number  placing  a  call  in,  say  (0,  t]  , 
with  t  <;  T  ,  has  the  Binomial  distribution  with  parameters  n  , 
p(t/T)  .  Letting  this  random  variable  be  N.  ,  we  note  that 
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E(Nt)  =  At 

where 

A  =  np/T  , 


and  E  denotes  expectation.  How  let  n  approach  infinity  and  p 
approach  zero  such  that  A  ,  the  mean  number  of  calls  per  unit  time, 
remains  constant.  The  limiting  process  is  a  Poisson  process.  The 
distribution  of  N.  approaches  the  Poisson,  mean  At  ,  and  the  times 
between  successive  calls  are  independent  and  have  the  exponential  dis- 
tribution. One  can  imagine  that  a  large  population  of  potential  callers, 
each  one  of  whom  has  a  very  small  probability  of  calling  during  (0,  T]  , 
and  all  of  whom  act  independently,  would  approximate  this  limiting  case 
quite  well.  The  Poisson  process  is  commonly  referred  to  as  the  model 
for  complete  randomness  in  a  point  process. 

Another  description  of  the  Poisson  process  in  terms  of  N,  will 
help  to  illustrate  the  model.  Note  that  w.  ,  the  number  of  events 
during  (0,  t]  ,  has  a  Poisson  distribution,  as  this  is  the  limit 
distribution  for  a  Binomial  with  increasing  n  ,  decreasing  p  and  con- 
stant mean  np  .  Define  N.  ,+,   to  be  the  number  of  events  occurring 
in  (t,t+h]  ,  so  that 

N+  4,,.  -  N.  _,.  -  U,  . 
t,t+h    t+h    t 

The  following  assumptions  also  characterize  the  Poisson  process. 

(a)  N.  .+,   has  the  Poisson  distribution  with  parameter  Ah  ,  for 
all  positive  t  and  h  . 
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(b)  N.  ...   is  statistically  independent  of  the  number  and  position 
of  the  events  in  (0,  t]  ,  for  all  positive  t  and  h  . 
Some  immediate  consequences  of  (a)  and  (b)  are: 

(1)  Prob  (Nt  t+h  =  n}  =  Prob  {Nh  =  n}  =  e~Xh(Ah)n/n!  regardless 
of  t  ,  so  that  the  process  is  stationary  in  time. 

(2)  Events  occur  singly  since  Prob{N.  ^  2}  ~  (Ah)  /2  near  h  =  0  . 

(3)  In  a  short  interval  of  length  h  ,  the  probability  of  one  event 
occurring  is  approximately  Ah,  and  of  no  events,  1  -  Ah  . 

(4)  The  process1  history  prior  to  time  t  and  its  evolution  after 
t  are  statistically  independent. 

From  (1)  it  follows  that  the  first  event  occurs  after  h  with 

probability  e~   ,  so  that  X,  ,  the  time  to  the  first  event,  has  the 

exponential  distribution.  This  reasoning  applies  to  successive  interevent 

times  X   by  conditioning  on  the  value  of  t  ,  ,  so  all  interevent  times 
n  J  3  n-1 

have  the  exponential  distribution.  The  independence  of  these  times  also 
follows  from  (4). 

To  see  that  our  definition  is  consistent  with  assumptions  (a)  and  (b), 
we  note  that 

!Nt  <  n}  =  {tn  >  t}  . 

This  equivalence  relates  the  counting  variable  N.   to  the  epochs  t   and 

implicitly  to  the  interevent  times  X   since  t  =  X,  +  ..  +  X  .  For 
v  J  n        n    1        n 

the  Poisson  process,  the  X   are  iid  with  the  exponential  distribution, 
so  t   has  the  Gamma  distribution,  yielding 
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Prob{N+  :  n}  =  Prob{t„   t}    J  e~At(At)k/k! 


n-1 
k=0 


It  follows  that  N.  has  the  Poisson  distribution,  with  mean  At  .  The 
memoryless  property  of  the  exponential  distribution  is  used  to  show  (a) 
and  (b)  hold  [Ross,  1972,  p.  118]. 

The  Poisson  process  enjoys  many  properties  not  shared  by  other  point 
processes.  These  properties  account  in  part  for  the  wide  use  of  the 
Poisson  process  in  modelling  problems.  The  counting  variable  N,  .  .  , 
the  number  of  events  in  (t,t+h]  ,  is  of  primary  interest.  It  has  the 
Poisson  distribution  with  mean  Ah  and  becomes  nearly  normal  as  Ah 
becomes  large.  By  using  a  mean  of  Ah  and  standard  deviation  of  /Ah  , 
the  normal  approximation  is  quite  reasonable  when  Ah  ^  20  .  The 
counting  variables  corresponding  to  any  set  of  disjoint  intervals  are 
statistically  independent.  This  is  particularly  useful  when  accounting 
for  the  effects  of  events  in  different  time  periods,  as  the  joint 
distribution  of  counts  is  easily  specified. 

Any  (non-random)  number  of  independent  Poisson  processes  can  be  super- 
imposed on  a  common  time  axis.  The  result  is  a  Poisson  process  whose  rate  is 
the  sum  of  the  rates  of  the  contributing  processes.  This  follows  directly 
from  the  reproductive  property  of  the  Poisson  distribution  under  summation. 
It  allows  one  to  combine  several  processes  without  sacrificing  any  proper- 
ties.  In  fact,  when  many  independent  non-degenerate  point  processes,  not 
necessarily  Poisson,  are  combined,  the  pooled  process  usually  tends  to  a 
Poisson  process.  This  fact  also  accounts  for  the  general  acceptance  of  the 
Poisson  process  when  modelling  the  output  of  several  independent  processes 
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of  unknown  origin.  See  [Cox  &  Smith,  1953,  1954]  or  [Khintchine,  1960] 
for  a  general  discussion  of  the  superposition  of  many  processes  and  of 
its  approach  to  tne  Poisson  process. 

Another  operation  under  which  the  Poisson  process  is  invariant  is 
the  random  deletion  of  events.   Suppose  that  every  event  is  subject  to 
deletion  with  probability  1  -  p  ,  independent  of  other  deletions. 
The  remaining  events  form  a  Poisson  process  with  rate  pA  ,  where  X 
was  the  original  rate.  This  result  follows  by  either  showing  that  the 
new  interevent  times  are  geometrically  compounded  exponentials  and  hence  are 
still  exponential,  or  by  showing  properties  (a)  and  (b)  hold  by  simple 
conditional  probability  arguments  applied  to  the  counting  variable.  This 
situation  would  arise  if,  for  example,  the  original  process  could  not  be 
monitored  directly,  such  as  suicides  from  a  bridge,  or  if  one  only 
recorded  certain  events,  such  as  certain  colored  cars  in  a  traffic 
stream;  the  recorded  events  would  remain  a  Poisson  process.  As  in  the 
case  of  superposition,  stationary  point  processes  which  are  thinned 
this  way  become  indistinguishable  from  Poisson  Processes.  By  normalizing 
so  the  mean  interevent  time  is  one,  the  characteristic  function  for 
interevent  times  of  the  thinned  process  tends  to  1/(1-0)  ,  continous 
at  the  origin,  as  the  probability  p  of  event  observation  tends  to 
zero.  This  corresponds  to  the  exponential  probability  distribution  with 
unit  mean.  Hence,  series  of  events  with  substantial  random  loss  should 

follow  a  Poisson  process  approximately. 

th 

The  epoch  t   at  which  the  n    event  occurs  is  also  of  interest. 

As  noted  earlier,  t   has  a  gamma  distribution  with  mean  n/A  and 
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variance  n/A  .  This  is  also  equivalent  to  the  fact  that  2At   has 

tne  chi-squared  distribution  with  2n  degrees  of  freedom,  so  standard 

tables  can  be  used  to  find  percentiles.  This  fact  also  provides  a 

confidence  limit  or  interval  for  A  based  on  observation  of  t  . 

When  n  becomes  large,  the  central  limit  theorem  shows  that 

(At  -  n)//n  has  the  standard  normal  distribution  approximately.  For 

improved  accuracy,  note  that  /At   is  nearly  normal  with  mean  /n-1/4 

and  variance  1/4. 

The  conditional  joint  distribution  of  (t-,,...,t  )  ,  given  the  event 

{N.  =  n}  ,  is  the  same  as  the  distribution  of  the  order  statistics  from 

n  independent  samples  of  a  uniform  distribution  over  [0,t]  .  This  means 

that  by  conditioning  on  the  number  of  events  in  an  interval,  the  event 

positions  in  the  interval  have  a  relatively  simple  probabilistic  law. 

Several  statistical  tests  for  deviations  from  the  Poisson  process  are 

based  on  this  [Epstein,  I960].  Ross  [1970,  p.  18]  uses  this  property 

to  show  that  in  a  queue  with  an  unlimited  supply  of  servers  and  Poisson 

arrivals,  the  number  of  busy  servers  at  time  t  is  distributed  Poisson 

with  mean  n(t) , 

t 
n(t)  =  A  Jo  (1-G(x))dx  , 

where  A  is  the  arrival  rate  and  G  is  the  distribution  of  service 
times.  Since  for  large  t  the  integral  becomes  E(S)  ,  the  mean 
service  time,  the  steady  state  solution  has  the  number  of  busy  servers 
distributed  Poisson  with  mean  \E(S)  .  This  model  may  arise  in  the 


design  of  service  systems  where  the  number  of  servers  is  to  be 
determined.  The  unlimited  server  case  suggests  the  actual  number 
needed  for  a  low  delay  response.  It  may  also  arise  where  service  is 
unlimited,  such  as  for  modelling  the  number  recovering  from  industrial 
accidents. 

The  local  variables  of  a  point  process  refer  to  those  variables 
measuring  the  process  at  the  time  t  .  For  the  Poisson  process,  these 
variables  have  simple  laws.  The  forward  recurrence  time,   Y,  ,  is 
the  time  from  t  forward  to  the  next  event.  Since  Y.  exceeds  h 
when  N.  .  .  =  0  ,  it  follows  that  Y.   has  the  exponential  distribution 
with  parameter  A  .  The  backward  recurrence  time,  Z.  ,  is  the  time 
since  the  last  event  prior  to  t  ,  or  t  if  none  has  occurred.  By  the 
same  reasoning,  Z   has  the  exponential  distribution  with  parameter 
A  but  is  truncated  at  t  .  From  property  (b),  Y.  and  Z,   are  inde- 
pendent, so  that  the  span,  S.  =  Y.  +  Z.   has,  for  t  »  A~  ,  the 
gamma  distribution  with  mean  2X   .  This  means  that  the  interevent 
time  which  extends  across  the  fixed  time  t  is  not  a  typical  interevent 
time,  as  it  has  twice  the  mean  of  usual  interevent  times  and  a  different 
distribution.  This  is  a  potential  source  of  confusion  and  error  in 
modelling  problems,  and  it  should  be  recognized  at  the  outset. 

A  Poisson  process  in  multiple  dimensions  is  called  spatial  in  [Feller 
1968,  p.  159].  He  gives  several  examples  of  this  case  with  data. 
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I V .  RENEWAL  PROCESSES 

When  a  Poisson  process  is  generalized  by  allowing  the  interevent 

times  to  have  a  distribution  distinct  from  the  exponential,  the  result 

is  a  renewal  process.  To  fix  the  notation,  we  use  t  ,  0  <,   t-,  ^  t?  <  . 

as  the  time  of  the  n    event,  or  renewal.  The  interevent  times  are 

X  =  t  -  t  i   ^  0  ,  which  are  taken  to  be  independent  random  variables, 
n    n    n-1     '  K 

It  is  sometimes  convenient  to  allow  X,  to  have  a  distribution  distinct 
from  the  others;  however,  for  now  we  assume  not,  so  that 


Prob{Xn  ^   x}  =  F(x)  ,  all  n  , 

and    Prob(X  >  x)  =  F(x)  . 
n        v  ' 

We  also  assume  that  F  is  a  non-lattice  distribution,  except  in  the 
section  on  numerical  methods  for  the  discrete  case.  The  second  moment  of 
F  is  assumed  to  be  finite,  with 

X"1  =  f  xdF(x) 


and    a  =  {  x  dF(x)  -  A" 

J   0 


In  subsequent  formulas,  the  product  Xo   ,  referred  to  as  the 
coefficient  of  variation  c  ,  appears  as  a  correction  factor  to  the 
exponential  case.  It  makes  a  convenient  and  reasonable  measure  of  the 
variability  of  a  renewal  process. 

As  before,  we  will  use  N.   as  the  number  of  events  in  (0,t]  or 

N.  =  max{n|  t  <;  t}  . 
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For  problems  involving  N.  ,+,  =  N.+,  -  N.  ,  we  can  refer  to  a  renewal 
process  whose  first  interevent  time  is  Y  ,  the  forward  recurrence 
time,  as  demonstrated  later. 

From  the  central  limit  theorem  applied  to  the  sums  t  ,  and  the 
equivalence 

{Nt  <  n)  =  Un  >  t}  , 


we  have  that,  as  t  -»  °°  , 


N  -  At  d 

+  Normal  (0,1) 
Aa/At 


This  means  that,  while  the  distribution  of  N.   usually  is  intractable, 
a  normal  approximation  for  large  t  exists,  and  standard  tables  can  be 
used  to  obtain  percentiles.  Note  that  the  coefficient  of  variation  appears 
as  a  correction  factor  to  the  standard  deviation  for  the  Poisson  case. 

A  quantity  of  consideraDle  interest  in  modelling  with  a  renewal 
process  is  the  renewal  function. 

M(t)  =  E(Nt)  . 

This  function,  always  non-negative  and  non-decreasing,  appears  in  most 

formulas  and  is  equivalent  in  information  content  to  the  distribution 

F  .  Both  functions  share  common  properties  of  differentiability,  so  that 

if  F  has  a  density,  f  ,  then  M  has  a  derivative,  m  ,  for  which 

t 
M(t)  =  J0  m(x)dx  , 
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where    m   ,     the  renewal    intensity,  may  involve  impulse  terms  if     f 
does. 

A  preliminary  limiting  result  for  the  renewal    function   is  that 

M(t)   =  At  +  o(t) 
as     t  ■*-  °°  , 

where  o(t)  denotes  a  quantity  tending  to  zero  when  divided  by  t  . 
The  only  renewal  process  with  a  common  interevent  distribution  for  which 
M(t)  =  At  is  the  Poisson  process.  Under  some  regularity  conditions 
(see  [Smith,  1962]),  the  analogous  result  for  the  renewal  intensity 
holds,  namely 

m(t)  =  A  +  o(l) 

as  t  -*■  °°  , 

where  o(l)  denotes  a  quantity  tending  to  zero.  Both  of  these  results 
are  intuitively  suggested  by  the  definitions. 

For  convenience,  we  will  assume  that  the  density  f  and  renewal 
intensity  m  exist;  however  the  methods  illustrated  do  not  depend  upon 
this.  There  are  two  lines  of  reasoning  which  can  be  identified  that 
work  well  when  modelling  with  a  renewal  process.  The  first  involves 
conditional  expectation  methods  and  the  second  an  interpretation  of  the 
renewal  intensity.  We  will  illustrate  both  methods,  as  they  often 
complement  each  other.  The  "renewal  equation"  on  which  some  treatments 
of  renewal  theory  are  based,  is  discussed  first  (see  [Feller,  1941]). 
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To  derive  an  equation  for  the  renewal  function,  M  ,  first  condition 
on  X,  ,  the  time  of  the  first  renewal.  We  have  that 

(l  +  M(t-x)   x  <:  t  , 
E(N.|X1  =  x)  ={ 

(    0       x  >  t  . 

Unconditioning  by  tne  distribution  of  X-,  ,  we  see  that 

t 

M(t)  =  E(N.)  =  F(t)  +  /  H(t-x)f(x)dx  (1) 

t  o 

and  by  differentiation 

t 

m(t)  =  f(t)  +  J  m(t-x)f(x)dx  .  (la) 

o 

We  will  refer  to  this  as  a  renewal  equation;  several  other  equations 
which  follow  appear  different  but  are  equivalent.   By  taking  the  Laplace 
transform  of  this  equation,  (transforms  being  denoted  by  an  asterisk*) 
we  have 


or  equivalently 


nms)  -  sT0ifsT] 


■Ms)  -  r^tsi 


where  for  any  g,  g(t)  =  0   for  t  <  0  ,  g*(s)  =  /  e~s  g(t)dt. 

This  equation  serves  to  show  that  the  renewal  function  and  the 
interevent  distribution  are  equivalent  in  information  content,  and  hence 
suggests  that  equations  involving  a  renewal  process  can  be  expressed  with 
the  renewal  function.  As  we  shall  have  occasion  to  use  results  for  the 
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case  where  X,  has  a  distribution  different  from  the  other  interevent 
times,  we  leave  as  an  exercise  for  the  interested  reader  to  show  that 


M 


}(t)   =  F^t)  +  J  F1(t-x)m(x)dx 


(lb) 


where  quantities  for  the  general  case  are  subscripted  by  1. 

The  conditional  expectation  approach  can  be  used  to  obtain  the 
distribution  of  forward  recurrence  time.  As  earlier,  we  take  Y   to  be 

the  time  elapsed  from  t  until  the  next  event,  or  t  -  t  for  n  =  l\L  +  1 

n  t 


Defining 


G(y,t)  =  Prob{/t  <;  y} 


we  have  by  conditioning  on  the  time  until  the  first  renewal 

G(y,t-x)  x  ss  t  , 

Prob  {Yt  £  y      :  x}  =<1  t  <  x  <,   t  +  y  , 

0  t  +  y  <  x  . 

After  unconditioning,  we  have 


G(y.t)  =  /  G(y,t-x)f(x)dx  +  F(t+y)  -  F(t)  , 

o 

which  after  some  manipulation  using  the  renewal  equation  (1)  and 
transforms  takes  the  form 

t+y 
G(y,t)  =  /   F(t+y-x)m(x)dx  . 
t 

To  illustrate  the  alternative  method  of  derivation  mentioned  above, 

observe  that  m(x)6x  is,  to  first  order  in  6x  ,  M(x+6x)  -  M(x)  ,  the 


(2) 


mean  of  N 


x+6x  .  As  6x  +  0  ,  this  mean  is  ProbiN^^  =  1}  +  o(5x)  . 
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In  essence,  the  probability  of  multiple  renewals  in  (x,x  +  6x]  is  of  order 
(6x)   ,  so  at  most  one  renewal  need  be  accounted  for  and  the  mean 
number  to  occur  is  asymptotically  the  probability  of  exactly  one  renewal 
in  (x,x+6x]  ,  as  6x  >  0  ,  Hence,  m(x)iSx  can  be  viewed  as  the  pro- 
bability of  a  renewal  "at  x".  Equation  (la)  can  now  be  reasoned  by  seeing 
that  m(t-x)6t  is  the  (conditional)  probability  of  some  renewal  occurring 
at  t  ,  given  the  first  renewal  occurs  at  about  x.  From  the  law  of 

total  probability,  and  by  passage  to  the  limit  in  the  choice  of  sub- 

t 
intervals  for  x  ,  we  have  that  /  m(t-x)f (x)dx-6t  is  the  probability 

0 

of  some  renewal  other  than  the  first  occurring  at  about  t  ,  to  first  order 
in  5t.  Now  as  f(t)6t  is  the  probability  of  the  first  renewal  occurring 
at  about  t  ,  to  first  order  in  6t  ,  we  have 

t 

m(t)St  -  (f(t)  +  {  m(t-x)f(x)dx)6t  +  o(6t)  ; 

J  o 

upon  dividing  by  6t  and  letting  6t  -*■   0  ,  equation  (la)  follows. 

Equation  (2)  is  obtained  in  a  similar  manner,  by  observing  that 
F(t+y-x)m(x)6x  ,  to  first  order,  is  the  probability  that  the  last  renewal 
prior  to  t+y  occurs  at  about  x  ,  and  that  this  event  occurs  at  some 
x  beyond  t  exactly  when  Y.  <;  y  .  The  case  in  which  no  density  for 
the  interevent  distribution  exists  is  handled  by  using  the  Lebesgue  - 
Stieltjes  integral  and  by  replacing  m(x)dx  by  dM(x)  . 

Obtaining  values  for  the  distribution  G(y,t)  depends  on  being  able 
to  obtain  the  renewal  intensity.  In  a  later  section,  some  approximations 
are  discussed.  Here,  the  steady  state  result  for  large  t  is  given: 
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y  _ 

G(y.t)  -  Fp(y)  =  A  {  F(x)dx  (3) 


as  t  -*■  °°  , 


where  F   is  called  the  equilibrium  distribution  corresponding  to  F  . 
This  result  can  be  obtained  from  equation  (2)  by  seeing  that,  for  large 
t  ,  m(x)  =   A  in  the  range  of  integration  (t,t+y]  .  Replacing 
m(x)  by  A  in  (2),  we  see  that,  as  t  ->  °°.  , 

t+y_ 

G(y,t)  =  j       F(t+y-x)Adx  =  F  (y)  . 
t 

A  rigorous  limiting  arguement  is  the  Key  Renewal  theorem  which  is  proven 

in  [Smith,  1958]. 

When  the  renewal  process  is  first  observed  at  a  time  t  from  its 

actual  origin,  the  first  interevent  time  is  Y,  .  For  many  practical 

situations,  one  must  assume  t  is  large  enough  to  use  steady  state 

results  for  determining  the  distribution  of  Y,  .   If  we  let  X,  have 

the  equilibrium  distribution  F   for  F  ,  while  all  the  remaining 

interevent  times  have  distribution  F  ,  then  we  have  that 

t 
Me(t)  =  Fe(t)  +  /0  Fe(t-x)m(x)dx  , 

where  M  (t)  is  the  mean  number  of  renewals  in  (0,t]  for  this  equilibrium 
renewal  process.  Taking  the  Laplace  transform  of  this  equation  gives 

M*(s)  =  F*(s)  +  F*(s)m*(s) 

and  since  *  ^ 

Fe(s)   ^(1  -  f  (s))  , 
e     s 
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*      A 
we  have  M  (s)  =  2  » 

e  '   s 


or 


Me(t)  =  At 


Thus  for  the  renewal  process  which  is  observed  in  steady  state,  the  renewal 
function  is  At  . 

The  backwards  recurrence  time  distribution  can  be  obtained  in  the 
same  manner  as  in  the  forward  recurrence  case.   In  particular,  the  same 
limiting  distribution  is  obtained; 

ProbiZ.  ^  z)  >  F  (z) 
t       e  ' 

as  t  ->  °°  , 

where  Z.   is  the  time  since  the  last  renewal  prior  to  t  .  Note  that  the 

exponential  distribution  is  its  own  equilibrium  distribution,  and  is  unique 

in  this  respect.  Finally,  we  note  that  our  previous  use  of  the  renewal 

intensity  immediately  yields  F(t-z)m(z)  as  the  density  of  Z   in  (0,t)  , 

with  probability  mass  F(t)  at  z  =  t  .   Integrating  this  density  gives 

t 
1  =  F(t)  +  /  F(t-z)m(z)dz  (4) 

0 

which  is  equivalent  to  the  renewal  equation  (1). 

The  forward  and  backward  recurrence  times  sum  to  the  span,  S.  , 
which  is  the  length  of  the  interevent  interval  intersecting  t  .  As  noted 
in  the  Poisson  process  case,  the  span  does  not  have  the  interevent  distri- 
bution F  .  The  sampling  of  interevent  intervals  is  biased  toward  the  longer 
intervals,  as  these  are  more  likely  to  intersect  the  fixed  time  t  .   If 
h(s,t)  is  the  density  of  S.  ,  then  for  s  <  t  ,  the  probability  h(s,t)  •  6s 
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is  the  probability  of  a  renewal  occurring  at  x  and  the  subsequent 
interevent  interval  being  of  length  s  ,  summed  over  x  in  [t-s,t]  , 

or  t 

h(s,t)  =  /   f(s)m(x)dx  =  f (s)[M(t)-M(t-s) J  . 
t-s 

As  t  tends  to  infinity,  h(s,t)  becomes  h(s)  ,  where 

h(s)  =  Asf(s)  , 

the  steady-state  density  of  span.  This  result  is  obtained  by  observing 

that  M(t)  ~  At  for  large  t  .  The  result  is  actually  an  application  of 

Blackwell's  theorem,  [Blackwell,  1948]  .  The  limiting  result  can  also  be 

obtained  by  arguing  that  in  a  set  of  n  random  intervals,  X,,...,X  , 

the  combined  expected  length  of  all  those  of  size  x  such  that  |x-s|  <  6_s_  , 

n  2 

is  approximately  n  •  s  •  f(s)  •  6s  which  normalized  by   £  X.  is  the 

i=l   ] 
probability  of  a  span  being  within  6s  of  s  .  As  n  becomes  large,  we 

have  the  limiting  result  derived  above  for  the  asymptotic  density  of  span 

by  the  strong  law  of  large  numbers. 

The  remaining  paragraphs  of  this  section  are  devoted  to  the  description 
of  examples  of  applications  involving  the  renewal  process  model. 

An  early  application  of  the  model  is  to  the  theory  of  counters,  as  in 
[Feller,  1948].  A  Geiger-Mul ler  counter  registers  the  emission  of  particles 
from  a  decaying  substance,  however  due  to  the  counter's  resolving  time,  not 
all  events  are   recorded.  The  evolution  of  recorded  events  can  be  modelled 
as  a  renewal  process  from  which  an  estimate  of  the  original  process  parameters 
can  be  made. 
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Models  for  a  reliability  problem  frequently  involve  a  renewal  process. 
See  [Gnedenko  et.  al . ,  1969]  for  a  treatment  of  reliability  theory  using 
renewal  theory.  The  interevent  intervals  correspond  to  machine  operating 
periods  and  the  renewal  epochs  to  machine  failures.  The  tail  distribution 
for  forward  recurrence  time,  ProbiY,  >  y}  ,  is  the  interval  reliability  as 
defined  in  [Barlow  and  Proschan,  1965,  p.  8].  A  model  for  the  single 
item  failure  and  repair  process  is  the  alternating  renewal  process.   In 
this  model,  the  process  is  in  one  of  two  states,  0  and  1  ,  and  alternates 
between  them.  Residence  times  in  state  0(1)  are  identically  distributed 
samples  from  a  distribution  F0(F-.)  >  ancl  al  1  residence  times  form  an 
independent  set.  The  epochs  at  which  the  process  enters  state  0  ,  or 
state  1  ,  form  two  imbedded  renewal  processes.  The  interevent  distri- 
bution for  either  imbedded  process  is  simply  the  convolution  of  FQ  with 
F-,  .  To  find  the  availability  function,  A(t)  ,  which  is  the  probability 
that  the  process  is  in,  say  state  1  at  time  t  ,  assuming  it  began  in 
state  1 ,  we  have  that 

t 
A(t)  =  F^t)  +  Jo  [F1(t-x)]m(x)dx  . 

The  reasoning  is  that  either  state  1  was  not  left  in  [0,t]  or,  after 
reentering  state  1  at  time  x  ,  the  process  remained  in  state  1  beyond 
time  t  ,  for  some  time  x  in  [0,t]  .  Some  bounds  in  special  cases 
are  given  in  [Butterworth ,  1963];  exact  solutions  are  usually  not  available. 
Approximations  can  be  found  by  using  one  of  the  approximations  to  m(x) 
reported  below.  The  steady  state  availability  follows  immediately  by 
replacing  m(x)  by  A  ,  its  relevant  value  in  the  integrand  when  t  is 
large,  to  obtain  y-. / (p-,  +  Mn)  in  the  limit.  Here  \i.      is  the  mean 
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residence  time  in  state  i.  The  interval  availability,  A(t,t+h)  defined 

as  the  probability  of  continued  residence  in  state  i  over  the  interval 

[t,t+hj  for  an  alternating  renewal  process,  has  an  analogous  derivation 

leading  to 

t 
A(t,t+h)  =  F, (t+h)  +  {  [Fn(t+h-x)]m(x)dx  .  (5) 

I  o    I 

The  steady  state  value  in  this  case  is  just  [u,/(ui  +  uQ)]  [F-,  (h)]  , 
which  is  suggested  by  the  form  of  the  expression. 

The  study  of  replacement  policies  for  failing  equipment  may  involve 
the  renewal  model .  Several  bounds  on  the  renewal  functions  are  deter- 
mined in  [Barlow  and  Proschan,  1964]  (and  by  alternate  methods  in  later 
sections  here)  by  comparison  of  certain  replacement  policies.   In 
particular,  for  a  block  replacement  policy  in  which  items  are  replaced 
every     b  units  of  time,  and  at  failure,  the  expected  number  of  failures 
between  preventive  replacements  is  M(b)  .  Accurate  estimates  of  this  may 
require  knowledge  of  the  renewal  function  for  relatively  short  times 
b.  The  approximations  given  below  are  particularly  directed  toward  the 
behavior  of  the  renewal  function  near  the  origin. 

Bartholomew  [1959,  1963  ]  applies  the  renewal  model  to  problems  in 
labor  turnover,  especially  in  circumstances  where  steady  state-results 
are  not  of  interest.  For  example,  the  renewal  intensity  is  interpreted 
as  the  labor  turnover  rate  for  new  organizations,  and  the  interevent 
distribution  as  a  means  of  measuring  stability.  An  approximation  to  the 
renewal  intensity  for  its  transition  from  the  origin  to  its  steady  state 
value  is  necessary  for  predicting  the  turnover  effect.  Actual  steady  state 
doesn't  occur  within  the  period  of  model  applicability. 
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The  queueing  theory  literature  relies  on  the  renewal  process  model 
to  describe  the  input  process  to  many  service  systems.  Marshall  and 
Wolff  [1971]  use  bounds  on  the  renewal  function  to  get  bounds  on  the 
difference  between  customer  average  and  time  average  measurements  for 
some  queues.  A  related  field,  Inventory  Theory,  has  problems  involving 
a  renewal  process  in  the  evaluation  of  the  s  -  S  policy  (see  [Arrow, 
Karlin  &  Scarf,  1958]). 

The  following  sections  will  give  some  bounds  and  approximations  with 
practical  significance,  and  briefly  review  some  related  literature. 

V.  THE  RENEWAL  FUNCTION 

As  we  stated  in  section  IV,  a  quantity  which  continually  arises  in 
renewal  theory,  and  in  the  analysis  of  many  models  in  which  renewal 
processes  are  imbedded,  is  the  renewal  function,  the  expected  number  of 
renewals  which  occur  in  the  interval  (0,t],  considered  as  a 
function  of  t  .  We  defined  this  to  be  M(t)  for  an  ordinary  renewal 
process,  and  showed  that 

t 

M(t)  =  F(t)  +  /  M(t-x)f(x)dx  ,     t  ^  0  .         (1) 

o 

Another  important  representation  of  M(t)  is  obtained  by  noting 
that  the  counting  random  variable  (N.  +1)  is  a  stopping  time  in  the  sense 
of  Wald,  and  by  looking  at  the  time  to  the  first  event  after  t  ,  that  is 
t  +  Y.  ,  after  taking  expectations 

M(t)  =  At  +  Ay(t)  -  1  ,     t  &  0 ■  ;  (6) 
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where  y(t)  =  E[Yt]  . 

a)  A£proxjmaticms_  for  the_Renewal_Fu_nctio_n 
For  all  but  a  few  simple  processes  it  is  difficult  to  solve 
equation  (1)  for  M(t).  It  is  equally  difficult  to  determine  M(t)  in 
other  ways,  such  as  from  (6),  but  since  the  renewal  function  appears  in 
so  many  expressions  for  other  quantities  of  interest,  it  is  important  that 
we  can  at  least  find  approximations  to  M(t)  .  Early  results  emphasized 
asymptotic  theory  for  large  t  .  Clearly  M(t)  ^  0  and  it  is  non- 
decreasing.  The  simplest  non-trivial  result,  often  called  the  elementary 
renewal  theorem  (see,  for  example  [Cox,  1962])  is  that  for  any  renewal 
process  with  finite  mean  interevent  time  1/A  , 

t   M(t) 

lim   )  •'  ■>  A  , 

t  -*•  °° 

so  that  A  can  be  interpreted  as  the  "rate"  of  the  process. 

Clearly  M(t)  increases  "on  the  average"  at  rate  \    .  The  next 
result,  which  can  be  stated  in  a  number  of  forms,  shows  that  M(t) 
approaches  a  linear  function  with  slope  A  (recall  we  are  assuming  that 
F  is  non-lattice.  When  F  is  lattice,  consider  M(t)  to  be  defined 
only  at  multiples  of  the  lattice;  then  M(t)  is  asymptotically  linear). 
This  well-known  result  is  known  as  Blackwell's  Theorem,  and  it  states 
that  for  any  h  >  0  ,  with  the  restrictions  just  discussed, 

lim  [M(t+h)  -  M(t)]  f  Ah  . 

t  -*  °° 

A  similar  result  is  found  in  [Smith,  1957]  .  Recall  that  a   is  the 

2    2  2 
variance  of  the  inter-event  times,  and  c  =  A  o  ,  the  coefficient 

-  22  - 


of  variation  squared.  Then  if  either  F  is  non-lattice,  or  if  t 
is  restricted  to  points  on  the  lattice,  then 

2  , 
M(t)  =  Xt  +  C--J  +  °(])  • 

Let  1 

A^t)  =  At  +  £^I  .  (7) 

In  this  paper  A. (t)  will  be  used  to  denote  the  ith  approximation 
to  M(t)  .  Equation  (7)  can  be  used  as  an  approximate  formula  for  M(t)  . 
It  can  be  quite  accurate  for  "large"  t  ,  where  "large"  depends  on  the 
distribution  F  as  well  as  its  mean.  For  many  distributions,  for  t  =  4 
or  5  multiple  of  tne  mean,  (7)  is  quite  accurate.  However  if  F  is  highly 
skewed,  convergence  to  the  linear  function  can  be  very  slow. 

Equation  (7)  has  certain  advantages  and  disadvantages.  First,  it 
is  linear,  easy  to  compute,  and  one  needs  knowledge  of  only  the  mean  and 
variance  of  the  interevent  time.  To  offset  these  advantages,  it  can  be 
quite  poor  for  approximating  M(t)  for  small  t  ,  and  it  is  not  known 

whetner  tne  true  M(t)  lies  above  as  below  A, (t)  for  any  given  t  . 

2-l 

Let  us  write  A,(t)  as  At  +  k  ,  where  k  =  ——■  .     Now  substitute 

this  approximation  to  M(t)  into  the  right  hand  side  of  the  renewal 
equation  (1)  .  After  a  little  algebra  the  expression  simplifies  to 
Ap(t)  ,  where 

A2(t)  =  At  +  (Hk)F(t)  -  Fe(t)  .  (8) 

This  gives  us  a  second  approximation  for  M(t)  .  It  is  easy  to  compute 
and  has  the  correct  asymptotic  behavior;  however,  we  have  lost  the  linear 
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property.  We  are  assured  that  A?(t)  is  non-negative.  This  last  property 

t  _ 
is  seen  by  noting  that  (1+k)  ;>  0  ,  and  F  (t)  =  A  j  F(u)du  ^  At  . 

*-  0 

From  equation  (1)  it  is  easy  to  derive  the  following  expression  for 
M(t)  eitner  directly  or  by  transform  methods, 

M(t)  =  At  +  /Q  Fe(t-u)m(u)du  -  Fe(t)  .  (9) 

Equation  (9)  is  useful  for  determining  another  approximation  to  M(t)  . 
Let  us  replace  m(u)  on  the  right  hand  side  of  (9)  by  A  ,  an  exact 
expression  for  the  case  of  a  Poisson  Process,  but  in  general  an  approxi- 
mation. Then 

t  _ 

A,(t)  =  At  +  A  /  F  (u)du  -  F  (t)  (10) 

3  J  o  e        e  x/ 

is  an  approximation  to  M(t)  .  Like  A~(t)  it  is  no  longer  linear,  and 
it  requires  knowledge  of  F  and  not  just  its  moments.  It  is  easy  to  show 

that 

c2  i 

lim  (A2(t)-At)  =  ■•--   , 

t  -*■   °° 

so  that  tne  approximation  is  asymptotically  correct. 

The  next  approximation  is  derived  from  an  approximation  to  the 

renewal  intensity  m(t)  derived  in  [Bartholomew,  1963].  Although  it  does 

not  appear  that  the  integral  of  the  approximation  was  intended  to  be  used 

as  an  approximation  to  M  ,  it  appears  to  give  wery   good  results  for 

small  t  ,  although  it  may  be  poor  for  large  t  .   In  the  notation  of 

this  paper,  Bartholomew  shows  that 

2 
m(t)  b  f(t)  +  A  ttfa    , 
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where  f  is  the  density  of  F  .  After  a  little  algebra,  the  integral  of 
his  approximation  (which  we  call  A*(t))  can  be  written 


A-(t)  =  U  +  F(t)  -  A  J 


1  - 


F(u)' 
FefuJ 


du 


(ID 


For  many  distributions  it  is  computationally  more  difficult  to  use 
(11)  than  (10)  or  (8)  because  of  the  quotient  in  the  integral.  Note 

that  the  denominator  F   is  itself  an  integral  of  F  . 

e  a 


Let 


K(t)  =  / 


1  - 


F(u) 


2 


F  Cu) 
e 


=  J  Fe(u) 


du 
2 


TeTQT 


~  •  du  . 


Let 


Then 


A(t)  =  F  It)   -   F(t) 


t  _ 


F(u) 


Now 


K(t)  =  JQ  F(u)du  +  Jo  A(u)  p-^}  •  du  . 
if  F  (t)  >  F(t)  for  some  t  ,  then  A(t)p|-K  <  A(t) 


Similarly  if  F  (t)  <  F(t)  the  same  inequality  holds.  Note  also  that 


since  A(t)  =  F(t)  -  Fe(t)  , 


J  A(u)du  =  !^ 


Thus  unless  F  (t)  -  F(t)  V  t  ^  0  ,  then 


K(») 


3-c' 
2A" 


-  25  - 


From  (11 )  we  see  that 


A3(t)  =  At  +  F(t)  -  AK(t)  , 


and  so  if  F  (t)  f  F(t)  for  some  t  ,  then 

2  . 
11m  [A3(t)-At]  >  ™-  -   1  . 

t  ->  °° 

Tne  approximations  given  so  far  are  not  restricted  to  any  particular 
family  of  distributions.  Papers  have  appeared  where  F  is  restricted  to 
a  particular  class.  Leadbetter  (1963]  proves  the  following  result.  Suppose 
we  can  write  F  as  an  absolutely  convergent  power  series 

r/4.\    v  (-1 )      .(tin      n 
FU)  ~~   J]  iVl+mn")  V   'm  >   °  ' 


Define 


Di  =  ci 


D2  =  C2  -  D,C, 


n-1 

d    =  c    -    y    D.C 

n    n   >.,  j  n-j 


Then 


00  /  i  ^n"l 
M(t)  =  J,  rlHvn]  Dnt   •  ■  >  °  •  (W) 


where  the  series  converges  absolutely  for  all  t  ^  0  . 

Applications  of  (12)  to  the  case  when  F  has  a  Weibull  distribution 

are  given  in  [Leadbetter,  1963]  and  [Smith  and  Leadbetter,  1963]  .  In  the 
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cases  treated,  only  the  first  four  terms  of  (12)  were  required.   It  should 
be  mentioned  that  they  treat  higher  moments  of  N,   and  in  fact  plot  the 
variance  of  N.  for  various  values  of  the  parameter  in  the  Weibull  dis- 
tribution. We  do  not  pursue  higher  moments  in  this  paper. 

Further  work  on  renewal  processes  for  which  F  has  a  Weibull 
distribution  is  to  be  found  in  [Lomnicki  ,  1966].  He  discusses  both  the 
distribution  of  the  number  of  renewals  and  the  renewal  function.  We  quote 
only  his  approximations  to  the  renewal  function. 

The  basic  idea  in  the  work  of  Lomnicki  is  to  write  the  renewal  function 
as  a  power  series  whose  terms  are  easy  to  calculate  from  available  tabu- 
lated functions,  and  where  convergence  is  rapid  so  that  only  a  few  terms 
need  to  be  calculated. 

Assume  that  the  underlying  distribution  F  is  given  by 

-t3 
F(t)  +  1  -  e    ,    ts:0,    3>0. 

Define  £L(t)  =  e_t  £  J-,  ,     k  =  1,2,3,  ...   . 

r=k 

The  functions  D  are  tabulated  as  the  tail  distribution  of  a  Poisson 
random  variable  with  mean  t  .  Now  define 

vfr\  _  r(3r-l)    r  o  i   2 

where  3  is  the  Weibull  parameter.  Next  define 

bQ(s)  =  y(s)  ,  s  =  0,  1,  2,  ... 

s-1 

bk+1(s)  =  I     bk(rb(s-r)  ,   k  =  0,  1, , 

r=k 

s  =  k+1,  k+2,  ...   . 
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Next  define 


ak(s) 


s 

I  (-1) 

p=k 


p+k 


bk(p) 

y(pT 


k  =  0,  1,  2, 
s  =  k,  k+1 , 


Finally  let 


ak(k)  =  ak(k)  , 


s        s-1 
(s)  =  I     a  (s)  -  y  af(s)  ,  s  =  k  +  1,  k+2, 
r=k        r=k 


Then  the  renewal  function 


M(t)      Ds(t3)   a.(s 


s 
s=l  °    k=l 

The  interested  reader  should  consult  the  Lomnicki  paper  for  details  such  as 
convergence  properties.  The  paper  shows  numerical  calculations  for  a  number 
of  interesting  functions  in  a  Wei  bull  renewal  process. 

b)  Bounds  for  the  Renewal  Function 

Instead  of  approximations  of  the  type  presented  in  a)  it  is  often 
desirable  to  nave  upper  and  lower  bounds  on  the  renewal  function,  and  a 
number  of  such  bounds  are  presented  here.   It  is  usually  easier  to  find 
lower  bounds.  Host  upper  bounds  apply  only  with  restricted  classes  of 
distributions,  although  some  general  upper  bounds  are  given. 

The  simplest  non-trivial  lower  bound  on  M(t)  is  obtained  from 
equation  (6).  By  noting  that  y(t)  is  the  expectation  of  a  non-negative 
random  variable  we  have  our  first  lower  bound  L-,(t)  ,  where 


L^t)  =  At  -  1 


(13) 


L.  will  be  used  for  the  ith  lower  bound,  and  U.  for  the  ith  upper  bound 
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A  lower  bound  was  derived  in  [Barlow  and  Proschan,  1965,  p.  68] 
in  an  interesting  way,  by  comparing  block  and  age  replacement  policies  in 
a  reliability  theory  context.  The  same  bound,  however,  is  easily  obtain- 
able from  equation  (9)  directly.  Since  F   is  a  monotone  non-increasing 
function,  equation  (9)  gives 

M(t)  £  At  +  Ffi(t)M(t)  -  Fe(t)  .  (14) 

Thus 

M(t)  £  L2(t)  =  Fn(t)     -  1  .  (14a) 

e^  ; 

As  was  pointed  out  earlier,  F  (t)  <,   At  ,  and  so  L?(t)  -•  0  .   It  is 
clear  also  that  since  F  (t)  $1    ,  that  L?(t)  ^  L,(t)  . 

Equation  (1)  can  be  used  to  find  a  simple  upper  bound  on  M(t)  .  Since 
F  is  monotone  non-decreasing  it  is  easy  to  see  from  (1)  that 

M(t)  <:  lUt)  =  U&   .  (15) 

1      F(t) 

In  general  this  bound  is  very  poor  and  is  of  little  practical  use.  The 
search  for  a  good  upper  bound  for  a  general  process  is  more  difficult. 
Stone  [1972]  in  a  short  interesting  paper  derives  the  following  very  useful 
upper  bound  which  has  the  nice  properties  of  being  linear  in  t  and 
requiring  only  the  first  two  moments  of  the  underlying  distribution; 

M(t)  <  U2(t)  -  At  +  3(l+c2)  .  (16) 

(Stone  in  this  and  a  number  of  other  papers  deals  with  renewal  processes 
in  which  the  random  variables  can  be  negative.  Since  our  motivation  is 
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taken  from  reliability,  inventory,  and  queuing  theory  where  the  random 
variables  are  invariably  non-negative  we  do  not  pursue  this  two-sided 
general ization. ) 

If  one  is  willing  to  restrict  the  class  of  distribution  functions, 
other  interesting  bounds  can  be  obtained.  Many  classes  of  distributions 
are   described  in  [Marshall  and  Proschan,  1970],  a  paper  which  studies 
the  relationship  between  various  classes. 

Let  us  assume  that  the  interevent  distribution  has  the  property 

/"tiujdu  ^  1/A  ,   all  t  ^  0  .  (17) 

t  F(t) 

The  left  hand  side  of  this  inequality  can  be  thought  to  represent  the 
mean  residual  lifetime  of  a  component  which  has  been  in  service  a  time 
t  and  which  has  a  random  lifetime  distributed  as  F  with  unconditional 
mean  1/A  . 

In  reliability  theory  terminology  we  say  that  F  belongs  to  the 
class  of  NBUE  distributions  (New-Better-than-Used-in-Expectation) ,  since 
clearly  the  inequality  says  that  a  used  component  can  never  have  a  mean 
residual  life  longer  than  the  mean  life  of  a  new  component.  Any  distri- 
bution for  which  the  inequality  in  (17)  is  reversed  is  called  UBNE  with 
the  obvious  connotation.  Equation  (17)  can  also  be  written  as 


Fe(t)  <.   F(t)  ,    all  t  ^  0  .  (17a) 


Equation   (1)   can  be  re-written  in  the  form 


t 
F(t)   =  J     F(t-u)m(u)du    .  (18) 


30 


Now  if  we  assume  that  F  is  NBUE(UBNE)  then  using  (17a)  and  (18)  in 
(9)  shows  that 

M(t)  ^  U  .  (19) 

( w 

In  this  paper  acronyms  (such  as  UBNE)  in  parentheses  should  be  read 
together  with  the  inequalities  in  parentheses.  Equations  (13)  and 
(19)  together  show  that  for  F  an  NBUE  distribution 

At  -  1  <;  M(t)  <:  At  ,  (20) 

and  we  have  very  simple  linear  bounds.  If  F  is  UBNE  (13)  can  be 
improved  to  K(t)  ^  At  . 

A  stronger  assumption  on  F  is  to  say  that  F  has  increasing 
failure  rate  ( I FR ) .  For  details  of  IFR  or  other  classes  the  reader 
should  consult  [Barlow  and  Proschan,  1965]  or[Marshall  and  Proschan, 
1972].  Here  we  shall  say  that  F  is  IFR  if  log  F(t)  is  concave  in 
t  .  For  a  renewal  process  with  F  restricted  to  this  class,  Barlow  and 
Proschan  [1965,  p.  71],  show  that 

M(t)  ^  At  pVt)   =  U4(t)  . 
e^ 

This  result  is  obtained  by  comparing  different  replacement  policies.  No 
direct  way  of  reproducing  this  result  has  been  found. 

The  linear  bounds  in  (2)  are  very  useful  in  many  applications.  Some 
are  discussed  later  in  this  paper.   If  each  bound  is  iterated  in  the  basic 
renewal  equation  (1)  ,  the  reader  will  find  that 
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L3(t)  =  At  -  Fe(t)  ^  M(t)  <;  At  -  Fe(t)  +  F(t)  =  U3(t) 


It  is  interesting  to  note  that  the  convex  combination  of  these  bounds 
which  has  the  correct  asymptotic  behavior  is  the  approximation  A«(t) 
in  equation  (8)  . 

Continued  iteration  of  either  bound  yields  a  sequence  of  bounds 
which  monotonically  converges  on  M(t)  from  below  or  above.  From  a 
computational  viewpoint  this  is  an  attractive  property.  The  question 
naturally  arises,  can  one  start  with  better  bounds  than  those  in  (20)? 
Of  course  the  upper  bound  in  (20)  holds  only  for  NBUE  distributions, 
but  Stone's  upper  bound  could  be  used  to  start  the  iterative  procedure. 

Marshall  [1973]  studied  the  problem  of  improving  the  starting 
bounds  in  this  iterative  procedure.   It  is  easy  to  show  that  for  any 
constant,  say  x  ,  if  the  function  At  +  x  is  used  to  start  an  iterative 
procedure  in  (1),  the  procedure  will  eventually  converge  on  M(t)  . 
Our  problem  is  to  find  two  numbers,  b   and  b   such  that 


At  +  b^   <;  M(t)  ^  At  +  ba  , 


where  Xt  +  b   is  as  large  as  possible,  and  At  +  b   is  as  small  as 
possible,  and  such  that  when  either  one  is  iterated  in  the  renewal 
equation  (1),  an  improved  bound  is  obtained  for  all  t  .  These  bounds 
are  called  the  "best"  linear  bounds. 

To  illustrate  the  method  we  use  the  lower  bound.  A  simple  iteration 
of  At  +  b.  in  (1)  shows  that  the  inequality 

t  By  "improved",  we  mean  a  bound  which  is  no  worse  than  the  one  in  question, 
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At  +  b£  <;  At  +  (l+b£)F(t)  -  Fe(t) 


must  hold  if  the  improvement  is  taking  place 
Thus 


FeU) 
*  F(t) 


bn  <;  Z**±i    -   1  , 


and  we  make  b  as  large  as  we  can.  Let  A  be  the  set  of  t  for  which 
F(t)  >  0  .  Thus  the  reader  should  see  that  the  largest  b  and  smallest 
b   are  given  respectively  by 

F  (t) 

b0  =  Inf  _^ -  1  , 

1       teA  F(t) 
and 

hi*) 

b  =  Sup  -  1  . 

U   teA  F(t) 


Since  Fg(0  )  =  F(0")  =  1  ,-1'sb.sO  .  Also  if  F  is  NBUE,  then 
F  (t)  <;  F(t)  ,  and  so  bu  =  0  .  In  this  case  (19)  is  indeed  the  "best" 
linear  upper  bound.  For  some  distributions  b   is  not  finite,  in  which 
case  M(t)  cannot  be  bounded  for  all  t  by  a  linear  function. 

c )  Examples  and  Calcula t i o ins . 

For  the  Poisson  process  (i.e.,  F  is  exponential)  the  solution  of 
(1)  is  M(t)  =  Xt  .  In  this  section  (1)  is  solved  for  a  number  of  renewal 
processes  and  compared  to  the  approximations  and  bounds  discussed  in  V 
a)  and  b) . 

The  first  case  considered  is  when  F  has  a  uniform  distribution  over 
(0,  2/A)  .  For  this  case  (1)  can  be  solved  for  successive  intervals  of 
length  2/X  ,  and  we  find 
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j-1  A^  -  n 

M(t)  p  y  i  (n  -  A2t)n  e2     -  1  , 
n=0  ' 

2(j-l)/A  ^  t  ^  2j/A, 

J  =  1,  2,  3,  ...   . 

Note  that  the  uniform  distribution  is  in  the  NBUE  class,  and  it  is  easy 
to  see  that  the  best  linear  bounds  are  given  by  (20)  .  Table  1  shows 
some  of  the  approximations  and  bounds  for  the  case  \   =  2  for  t  up  to 
10  mean  lifetimes. 

The  second  case  shown  is  when  F  has  a  gamma  distribution  with  mean 

2 
1  and  variance  -5  .  As  an  example  with  C  >  1  a  third  case  is  shown 

with  F  a  hyperexponential  distribution  with  mean  1  and  variance  2.5  . 

Calculations  for  the  gamma  and  exponential  cases  are  shown  in  tables  2 

and  3  respectively.  The  gamma  distribution  is  NBUE  with  k  =  -  .25  , 

b  =  -  .5  and  b  =  0  .  The  hyperexponential  distribution  is  UBNE  with 

k  =  .75  ,  b„  =  0  and  b  =  1.5  . 

I  u 
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V I •  THE  RENEWAL  INTENSITY  AND  ITS  CONVOLUTIONS. 

The  previous  section  gives  approximations  and  bounds  on  the  renewal 
function.  In  this  section,  some  approximations  for  the  renewal  intensity 
are  shown.  Also  some  approximations  to  convolutions  of  tail  distribution 
functions  with  the  renewal  intensity  are   given  and  some  applications  for 
these  are  discussed. 

The  renewal  intensity  function  m(t)  gives  the  rate  of  renewals 
occuring  at  the  instant  t  from  the  start  of  a  renewal  process.   If  many 
(say  N)  identical  renewal  processes  are  operating  simultaneously,  then 
N  •  m(t)  is  the  rate  of  renewals  for  the  composite  process.  This  fact 
is  at  the  basis  of  an  application  to  labour  turnover.  A  new  organization 
of  fixed  size  experiences  a  turnover  in  its  labour  force  which  has  been 
observed  to  fit  this  renewal  model.  See  [Bartholomew,  1959]  for  a  complete 
discussion  with  some  data. 

As  mentioned  in  an  earlier  section,  m(t)  approaches  under  suitable 
regularity  conditions,  the  rate  A  for  the  process,  as  t  tends  to 
infinity.  Our  approximation  is  meant  to  describe  the  transition  from 
M(0)  =  f(0)  to  the  final  value  of  A  .  Perhaps  the  best  approximation 
we  know  of  is  the  one  given  in  the  reference  [Bartholomew,  1963a]  , 
which  we  call  h,(t)  ,  where 

h  (t)  =  f(t)  i   .^M.   .  (21) 

1  J  F(u)du 

o 


As  seen  later  in  numerical  illustrations,  h,   is  an  excellent 
approximation  to  m  ,  particularly  for  highly  skew  interevent  distributions, 
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which  is  exactly  when  approximations  are  most  needed.  This  is  explained 
in  part  by  the  observations  given  by  Bartholomew  which  are  repeated  here. 
The  renewal  equation  (4)  derived  in  section  IV  is  essentially 


t 
1  =  F(t)  +  /  m(t-z)F(z)dz  , 


or  equalently 


t 
1  =  F(t)//  m(t-z)F(z)dz 


From  the  renewal  equation  (la)  we  also  have 

t 


m 


(t)  =  f(t)  +  /  m(t-u)f(u)du  .  (la) 


Upon  combining  these,  this  renewal  equation  may  be  written 

t  t 

m(t)  =  f(t)  +  F(t)  /  m(t-u)f(u)du//  m(t-u)F(u)du 


The  approximation  h,  then  results  by  replacing  the  ratio 

t  t 

/  m(t-u)f(u)du/J  m(t-u)F(u)du 

J  o  o 

by 

t       t  _ 

J  f(u)du/f  F(u)du  . 

J  o         o 

Since  m(t)  =  A  a  constant  for  the  exponential  case,  it  follows  that 
h-,  =  m  for  a  Poisson  process.  Further,  in  cases  where  F  is  highly 
skewed  so  that  m(t-u)  in  more  nearly  constant  in  u  where  f(u)  and 
F(u)  have  appreciable  positive  value,  the  factoring  and  cancellation 
of  m(t-u)  is  not  so  unreasonable  and  the  approximation  is  close. 
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Other  reasons  for  the  success  of  h,  are: 

(1)  ^(0)  ■  m(0) 

(2)  Lim  h^t)  =  Lim  m(t)  =  A 

t  "*■  °°       t  ■*■  °° 

(3)  h-j(O)  =  m'(0)  and  h^'(0)  =  m"(0) 

where  prime  denotes  differentiation  with  respect  to  the  single  variable  t. 

An  advantage  of  h-,  is  that  knowledge  of  the  interevent  distribution 

F  is  necessary  only  up  to  t;  that  is,  h,(t)  is  independent  of  the  tail 

of  F  beyond  t  .  This  can  be  of  considerable  advantage  when  F  is  skewed 

and  must  be  estimated  from  data.   In  particular,  methods  depending  on 

transforms  are  inherently  sensitive  to  the  tail  of  F  ,  a  function  which 

may  not  be  well  known  in  the  upper  regions,  particularly  if  data  is  not 

available.  The  only  possible  drawback  to  h,  is  the  appearance  of 
t  _  ' 

/  F(u)du  which  might  have  to  be  obtained  by  numerical  integration.  Even 

this  however,  should  not  be  difficult. 

Another  approximation  to  m(t)  to  consider  is  h«{t)  where 

h2(t)  =  f(t)  +  XF(t-)  .  (22) 

This  approximation  results  by  replacing  m(t-u)  by  A  in  the  right- 
hand  side  of  equation  (la)  .  Its  relative  lack  of  accuracy  when  com- 
pared to  h,  is  somewhat  compensated  for  by  the  simpler  functional 
form,  however  the  choice  of  which  one  to  use  is  best  left  to  those  making 
the  application.  While  h«  can  itself  be  substituted  in  the  righthand 
side  of  (la)  to  give  an  improved  approximation,  the  appearance  of 
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convolutions  quickly  limits  the  practable  usefulness  of  this  approach, 
particularly  in  view  of  the  availability  of  h, . 

Tables  4  &  5  shown  below  give  an  indication  of  the  numerical 
accuracy  involved  with  these  approximations. 


t 

m(t) 

h^t) 

h2(t) 

m(t') 

h^t) 

h2(t) 

0.0 

1.60 

1.60 

1.60 

4.51 

4.51 

4.51 

0.2 

1.51 

1.51 

1.36 

4.11 

4.11 

2.24 

0.4 

1.44 

1.44 

1.21 

3.75 

3.77 

1.40 

0.8 

1.32 

1.32 

1.04 

3.16 

3.25 

.99 

1.2 

1.23 

1.24 

.98 

2.69 

2.87 

.93 

1.6 

1.17 

1.18 

.95 

2.32 

2.60 

.93 

2.0 

1.12 

1.14 

.95 

2.04 

2.39 

.93 

2.5 

1.08 

1.11 

.95 

1.76 

2.18 

.94 

3.0 

1.05 

1.08 

.96 

1.56 

2.02 

.94 

5.0 

1.01 

1.03 

.98 

1.17 

1.62 

.95 

10. 

1.00 

1.00 

1.00 

1.01 

1.25 

.97 

A 

=   1.0   , 

c  =  1.58 

A    = 

=  1.0    ,     c 

=  3.54 

TABLE  4.  Calculations  using  a  Hyper-Exponential  interevent  distribution 
with  the  given  rate  A  and  coefficient  of  variation  c. 
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t 

m(t) 

hj(t) 

h2(t) 

0.0 

0. 

0. 

0. 

.01 

.039 

.039 

.039 

.02 

.077 

.077 

.078 

.04 

.148 

.148 

.151 

0.1 

.330 

.331 

.345 

0.2 

.551 

.556 

.598 

0.4 

.798 

.818 

.910 

1.0 

.982 

1.03 

1.14 

2.0 

1.00 

1.02 

1.05 

4.0 

1.00 

1.00 

1.00 

A  =  1.0   c  =  0.77 


TABLE  5.  Calculations  using  an  Erlang  interevent  distribution  with 
shape  and  scale  parameters  =  2.0  . 
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Another  renewal  process  quantity  to  be  approximated  is  a  con- 
volution of  the  form 

t 

c(t)  =  /  G(t-u)m(u)du 

o 

where  G(x)  =  1  -  G(x)  is  any  tail  distribution  function  and  m(u)  is 
the  renewal  intensity.  This  convolution  appears  in  application  formulas 
such  as  for  the  availability  in  an  alternating  renewal  process,  the 
mean  number  present  in  a  GI/G/  °°  queue,  and  the  distribution  of  the 
local  variables  forward  and  backward  recurrence  time.  Since  the  forward 
recurrence  time  has  the  widest  appeal  to  applications  generally,  we 
will  illustrate  both  proposed  approximations  in  the  context  of  an 
approximation  to  the  forward  recurrence  time  distribution. 

The  Key  Renewal  Theorem  assures  that,  since  the  interevent  dis- 
tribution F  is  assumed  to  be  non-lattice,  the  limiting  value  for  the 
convolution  in  question  is 

t 
lim  /  G(t-u)m(u)du  =  u-X  , 

t  -*-  co   ° 

where  °°  _ 

u  =  |  G(x)dx   (<co)  . 

o 

(Note:  u  is  the  mean  for  distribution  G  .) 

To  obtain  our  first  approximation  to  the  convolution  c(t)  ,  we 
can  simply  replace  m(u)  by  h,(u)  or  h2(u)  ,  and  use  numerical 
integration  to  obtain  values.  Since  h   is  generally  a  closer  approxima- 
tion, we  take 

t  _ 

c(t)  =  /  G(t-u)hn(u)du  (23) 

o        i 
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as  an  approximation  to  c(t)  .  As  we  seen  in  the  numerical  illustra- 
tions below,  this  can  be  a  very  close  approximation  when  used  for  the 
forward  recurrence  time  distribution. 

Another  approximation  to  c(t)  is  easily  obtained  as  follows.  By 
considering  G  to  be  the  distribution  of  X   in  the  renewal  process, 
and  M,  to  be  the  corresponding  renewal  function,  we  have  the  following 
version  of  equation  (lb) 

t  _ 

c(t)  =  /  G(t-u)m(u)du  =  G(t)  +  M(t)  -  M,(t)  .  (24) 

o  I 

By  replacing  M(t)  by  A«(t)  (see  equation  (8))  and  M,  (t)  by  an 
analogous  approximation  for  the  general  case,  we  have 


M(t)  ~  At  +  (Hk)F(t)  -  Fe(t) 


and 


M^t)  ~  At  -  Fe(t)  +  G(t)  +  [Hk-AyjF(t)  , 
which  upon  substitution  into  equation  (24)  yield 

c(t)  =  XuF(t)  (25) 

as  an  approximation  to     c(t)    . 

Applying  these  approximations   to  the  distribution  of  forward  recur- 
rence time,     Y     ,     we  have 

Prob{Y     >  y}  *  F(t+y)   +  J     F(t+y-u){f (u)   +  -JLM±  }du 

L  o  u 

f     F(x)dx 

0 
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as  the  first  approximation.  Values  are  obtained  by  using  numerical 
integration. 

To  apply  the  second  approximation  to  the  forward  recurrence  time 
distribution,  we  can  write 


where 


Prob{Y.  >  y}  =  F(t+y)  +  F(y)  J  G(t-u)m(u)du 


G(x)  =  F(x+y)/F(y)  , 


That  is,  G  is  the  conditional  distribution  of  residual  life  beyond  y 
In  this  case,  u  is  the  mean  residual  life,  or 


00  c/  _l  \     f  F(x)dx 

u  =  /  F-(-x^  dx  =  Jy,      -  . 

0  F(y)        F(y) 
Now,  substitution  of  equation  (24)  into  the  above  equation  yields 


Prob{Yt  >  y}  =  F(t+y)  +  Fe(y)F(t) 


The  following  table  illustrates  the  two  approximations  in  a 
particular  case. 


(26b) 
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0. 


0.2 


0.4 


0.6 


0.8 


1.0 


1.5 


2.0 


2.9 


Hyper-Exponential 
Interevent  Times 

2-Erlang 
Interevent  Times 

A  =  1.0  c  =  1.58 

A  =  1.0  c  =  .707 

t  =  0.5 

t  -  0.5 

1.0 
1.0 
1.0 

1.0 
1.0 
1.0 

0.765 
0.765 
0.804 

0.823 
0.827 
0.804 

0.600 
0.600 
0.663 

0.653 
0.657 
0.629 

0.483 
0.483 
0.558 

0.506 
0.509 
0.482 

0.398 
0.398 
0.479 

0.385 
0.388 
0.363 

0.335 
0.335 
0.419 

0.289 
0.291 
0.271 

0.236 
0.236 
0.314 

0.135 
0.135 
0.125 

0.179 
0.179 
0.247 

0.060 
0.060 
0.055 

0.119 
0.119 
0.167 

0.013 
0.013 
0.012 

TABLE  6.  Forward  Recurrence  Time  Distribution.  Entries  are  Prob{Y.  >  y}; 
1st  number  is  exact  value,  2nd  number  using  (26a),  3rd  using  (26b) 
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VII.  CALCULATIONS  IN  DISCRETE  RENEWAL  PROCESSES 

In  many  practical  applications  time  is  measured  in  discrete  units, 
and  the  assumption  of  continuous  distributions  is  an  approximation. 
Renewal  processes  find  application  in  inventory  theory  (for  example 
see  [Arrow,  Karlin  and  Scarf,  1958]  where  the  variable  t  is  not 
"time"  but  discrete  items  of  some  commodity.   In  this  section  we  present 
a  simple  computational  method  for  solving  renewal  equations  with  discrete 
distributions  and  show  some  computational  results. 

Let  X.  be  the  time  between  renewals  (i-1)  and  i  ,  and  let 


F  =  Prob{Xn  =  n}  ,    n  =  0,  1  ,  2,  ...  , 


be  the  probability  mass  function,  the  same  for  all  i  .  The  ideas  are 

first  illustrated  in  obtaining  the  solution  to  the  renewal  function, 

equation  (1).  In  the  discrete  case  let  M   be  the  expected  number 

n 
of  renewals  up  to  an  including  time  n  ,  and  F  =  J  =   ProblX  <;  n}  . 


Then  equation  (1)  can  be  written 


N  =  F  +  J  M   .F.  ,   n  =  0,  1,  2,  ...  . 
n    n   .^Q  n-j  j 

For  any  fixed  n  ,  define  the  (n+1)  x  (n+1)  lower  triangular  matrix 


F  = 


~fo 

0-  - 

-0 

1 

fl 

fo 

0 

f 

n 

f 

(1) 
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and  let  F  and  M  be  the  (n+1)  column  vectors 


and 


F  =  EF0,  Fr  F2,  ...  ,  Fn]  , 


M  =  [M0,  Mr  M2,  ...  ,  M  ] 


respectively.  Then  equation  (1)  can  be  re-written  as 


=  F  +  FM 


Note  that  F  is  a  non-negative  matrix,  every  row  sum  is  <  1  ,  and 
at  least  one  row  sums  to  strictly  less  than  1  (assuming  E[X]  >  0)  . 
Therefore  (I-F)~   exists,  where  I  is  the  identity  matrix,  and 


M  =  (I-F)"1 F 


The  inverse  matrix  of  (I-F)  has  the  structure 


G  = 


g0. 
% 


0 
90 


and  so  is  completely  determined  by  a  vector  g  =  [g  ,...,g  ]. 

For  small  n  G  is  probably  quite  easily  calculated  using  a 
general  matrix  inversion  routine.  The  primitive  operator  (jj  in  the 
APL  language  is  extremely  efficient.  But  for  large  n  it  is  advan- 
tageous to  make  use  of  the  special  structure  of  F  . 

Let  us  assume  that  (n+1)  is  a  power  of  2.   If  not  we  increase  it 
to  the  next  higher  power.  For  example  if  (n+1)  is  100  increase  it 
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to  128  =  2  .  Thus  let  us  assume  j  is  given  by  2J  =  (n+1),  and 

write  F,  for  the  (n+1)  x  (n+1)  matrix  F  ,  G.  for  the  matrix  G 
J  J 

If  we  partition  F .  into 

J 


F.   = 
J 


Fo-i     ° 


H,  ,  F.  , 


,   J   =   1,   2, 


with     FQ  =  F~     and     H-.  =  F,    ,     then 


'j+1 


G.  0 

J 


(Gj   Hj  Gj)Gj 


j   =  0,   1,   2,    ...    , 


with     GQ  -   (l/l-f0)    . 

Using  an  array  manipulation  language  such  as  APL  one  can  calculate 
large  dimensioned  matrices  G  in  only  a  few  iterations  with  no  matrix 
inversion  necessary. 

Fiqure  1  shows  the  renewal  function  for  X.  distribution  bino- 
mially  with  parameters  100  and  .2  .  Thus  E(X.)  =  20  ,  c  =  .2 
and  k  =  -.48  .  The  approximation  A,  (t)  =  .05t-.48  is  shown  for 
comparison.  A  four-instruction  program  written  in  APL  took  only  1.8 
seconds  on  an  IBM  370  model  158  to  solve  for  M  with  n  =  200  (10  mean 
times) . 

Having  determined  a  method  of  solving  (1)  we  have  really  solved 
a  whole  class  of  problems  in  renewal  theory.  In  place  of  F  let 
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FIGURE  1:  Renewal  Function__for_  a  Binomial  Process  with  Parameters, 
TOO  and  2. 
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us  write  a  vector  h  =  [hQ,  h-j ,  .  .  .  ,  h  ]  :  Now  consider  the  renewal 
equation  for  the  (n+1)  -  vector  r  , 

r  =  h  +  Fr  , 

where  F  is  the  same  matrix  as  before.  Then 

r  =  Gh  . 

By  carefully  choosing  h  a  number  of  problems  can  be  solved.  Let 

h  =  [fQ,  f i  ,  ...  ,  f  ]  .  Then  r  is  a  vector  of  renewal  intensities. 

Figure  2  shows  the  renewal  intensity  for  a  binomial  renewal  process 

with  parameters  100  and  .2  .  As  a  further  example  consider  the 

distribution  of  the  excess  random  variable.  Let  Y   be  the  time  until 

n 

the  first  renewal  after  time  n  ,  and  G(j,n)  =  Prob{Y  >  j}  , 
j  =  0,  1,  2,  ...   .  Define  for  a  fixed  j  G  =   [G(j,  0)  ,  G(j ,  1 )  ,  . . .  , 
G(j,  n)]  .  If  we  let  h  =  [F.,  r\+1 ,  F.+2,  ...  ,  F.+ J  ,  then  r  will 
give  us  G.  .  In  this  way  we  can  determine  for  a  fixed  j  how  the  excess 
distributions  are  converging  with  increasing  n  . 

Consider  an  alternating  renewal  process  as  described  in  Section  IV. 
Let  A   be  the  probability  the  system  is  operating  at  time  n  .  Let  the 
"up"  times  be  distributed  with  distribution  function  U.  =  P[U,  time  ^  j]  . 
Let  F.  be  the  probability  mass  function  of  the  sum  of  an  up  time  and  a  down 
time,  and  F  the  same  as  above.   If  h  is  taken  to  be  the  vector 
[Un,  U, ,  ...  ,  U  ]  then  r  will  give  the  availabilities  (assuming  the 
system  starts  in  the  "up"  state)  at  times  0,  1,  ...  ,  n  . 


-  51  - 


FIGURE  2:  Renewal  Intensity  for  a  Binomial  Process  with_ Parameters 
100  and  "72  . 
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VIII.  MISCELLANY. 

In  this  section  we  mention  a  few  topics  of  special  interest. 
These  topics  are  higher  moments  of  the  renewal  counting  variable,  some 
generalizations  of  the  renewal  model  and  superposition  of  several  renewal 
processes. 

In  some  applications,  higher  moments  of  the  renewal  counting  variable 
N,  may  be  of  interest.  In  particular,  the  variance  of  N.  may  be  desired 
Consider,  for  example,  a  situation  in  which  events  occur  at  times 
t.  =  i  •  a  +  e.  ,  where  e,  e~,  ...  ,  are  I  ID ,  normal  (0,  a)  ,  as  i 
runs  over  all  intergers.  This  may  be  the  arrival  pattern  of  persons 
with  appointments.  This  is  clearly  not  a  renewal  process  (unless  a  -   0)  , 
however,  when  a  is  large,  the  point  process  appears,  over  fixed  time 
periods,  similar  to  a  Poisson  process.  To  distinguish  this  point  process 
from  the  Poisson  process  or  perhaps  from  other  renewal  processes,  the 
variance  of  N.  can  be  investigated.  The  variance  of  N.   is  also  used 
in  the  analysis  of  the  pooled  output  of  several  renewal  processes,  a 
topic  mentioned  below. 

Higher  moments  of  N,  are  treated  in  [Cox,  1962,  pp.  59-60]  , 

where  the  Laplace  transform  methods  are  used.  The  r    semi-invariant 

moment  of  N   is  shown  to  be  asymptotic  to  A  t  +  V  +  0(1)  ,  (see 

also  [Smith,  1959])  where  A   is  a  function  of  the  first  r  moments 

of  the  interevent  distribution,  and  V   is  a  function  of  the  first 

r  +  1  moments.  In  particular,  the  variance  of  N   is  (\o)  At  +o(t) 

as  t  tends  to  °°  . 

Since  the  quantity  (Aa)   appears  as  a  factor  times  the  variance 
of  a  Poisson  Process  with  the  same  rate,  it  is  reasonable  in  practical 
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investigations  to  use  the  coefficient  of  variation,  Ao  ,  as  a  measure 
of  variation  in  the  counting  processes  Nt  .  Accordingly,  when 
Ao  <  1  (>  1)  ,  we  say  the  process  is  "under  variance"  ("over  variance"), 
the  boundary  or  bench  mark  case  being  the  Poisson  process  for  which 

Ao  =  1  . 

Some  bounds  on  higher  moments  are  provided  in  [Barlow  '6.   Proschan, 

1964],  where  it  is  proven  that 

E[ilt(Nt-l)  ...  (N  -k+1)]    fN j       (At)k 

— ki =  Eln  *  [*]    kv  ' 

when  we  have  that  F  is  NBUE  (UBNE)  respectively.   It  is  also  shown 
that  the  variance  can  be  bounded  by  the  renewal  function,  that  is, 

V(Nt)  <:  (*)  E(Nt)  =  M(t) 

when  we  have  that  F  has  increasing  (decreasing)  failure  rate  respectively, 

One  generalization  of  the  renewal  process  model  discussed  by  Smith 
[1958]  and  Cox  [1962]  is  the  cumulative  process.   In  this  model,  a  random 
amount  W.   is  contributed  at  each  renewal  epoch  t.  .  We  are  primarily 
concerned  with  the  net  contribution  made  by  time  t  .   If  Z.   represents 


this  net  contribution,  we  have 

N 


t 


Zt  =  I     Wi  ■ 
1   i=l  n 

where  the  empty  sum  is  given  value  zero.   If  the  W.  are  identically 

equal  to  one,  then  Z.  =  N,  .   If  W.   is  the  down  time  resulting  from  the 

i   replacement  of  a  failed  part,  Z.   is  the  down  time,  cumulative  to 

t  .  The  W.  may  assume  positive  or  negative  values,  if  for  example  W. 
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+■  h 

is  a  monetary  return  associated  with  the  i    event.  The  mean  and 
variance  of  Z  ,  when  the  W.  are  iid  ,  and  independent  of  N. 
are  given  in  [Cox,  1962,  pp.  94]  as 


E(Zt)  -  H(t)  ■  Mw 


V(Zt)  =  M(t)  •  o2m   +  V(Nt)„2 

2 

wnere  u.  and  a   are  the  mean  and  variance  respectively  of  the  W. 
WW  i 

distribution.  Asymptotic  normality  of  Z   is  also  shown,  with  the 
asymptotic  expansions  for  M(t)  and  V(N.)  replacing  the  exact 
values  in  the  above  expressions. 

In  [Barlow  &  Proschan,  1964],  the  generalization  in  which  the 
successive  times  between  events  may  have  different  distributions,  but 
must  have  a  common  mean,  is  discussed.   Interevent  times  continue  to  be 
an  independent  set.  The  bound 

M(t)  <  (>)  \t 

is  proven,  when  all  interevent  distributions  have  increasing  (decreasing) 
failure  rate  respectively.  These  hypotheses  can  be  relaxed  to  NBUE  (UBNE) 
respectively. 

The  superposition  of  several  renewal  processes  provides  a  point 
process  distinct  from  a  renewal  process,  (except  in  the  Poisson  Process 
case).  Such  processes  are  observed  when  two  or  more  arrival  streams  are 
pooled,  for  example,  or  when  failures  of  a  system  are   the  failures  of  any 
component,  and  each  component  fails  according  to  some  renewal  process. 
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The  asymptotic  result  for  many  independent  contributing  processes  is 
that,  over  intervals  of  fixed  length,  the  composite  process  exhibits 
the  properties  of  a  Poisson  process.  The  contributing  processes  may 
be  non-identical,  but  the  relative  contribution  to  the  total  rate  made 
by  each  contributing  process  must  tend  uniformily  to  zero.  See 
[Drenick,  1960]  for  a  complete  statement  of  tne  necessary  conditions. 
The  superposition  of  only  a  few  processes  is  treated  in  [Cox,  1962, 
pp.  71-77],  wnere  distributional  results  are  emphasized,  and  in  [Cox  and 
Lewis,  1966,  pp.  210-223]  where  a  statistical  analysis  is  given  as  well. 
See  also  [Cox  and  Smith,  1954]  for  a  complete  discussion. 
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